Apprentissage R

INTRODUCTION A LA STATISTIQUE AVEC R

I- Variables

---

Variables simples

x=2
y=6
2*x-y

[1] -2

Vecteurs

T=c(160, 180, 165, 190, 170, 140, 155, 167, 170, 169)
P=c(80, 75, 70, 85, 83, 87, 70, 89, 75, 85)
S=c("H","F", "F", "H", "H", "F", "F", "F", "H", "H")
C=c("bleue","noir","gris","bleue", "bleue", "gris", "bleue", "noir", "noir", "gris")

Matrices

df_mat=cbind(T,P)
df_mat

        T  P
 [1,] 160 80
 [2,] 180 75
 [3,] 165 70
 [4,] 190 85
 [5,] 170 83
 [6,] 140 87
 [7,] 155 70
 [8,] 167 89
 [9,] 170 75
[10,] 169 85

#NB: Les element d'une matrice sont de même type
cbind(df_mat, S, C)

      T     P    S   C      
 [1,] "160" "80" "H" "bleue"
 [2,] "180" "75" "F" "noir" 
 [3,] "165" "70" "F" "gris" 
 [4,] "190" "85" "H" "bleue"
 [5,] "170" "83" "H" "bleue"
 [6,] "140" "87" "F" "gris" 
 [7,] "155" "70" "F" "bleue"
 [8,] "167" "89" "F" "noir" 
 [9,] "170" "75" "H" "noir" 
[10,] "169" "85" "H" "gris" 

160

180

165

190

170

140

155

167

170

169

80

75

70

85

83

87

70

89

75

85

H

F

F

H

H

F

F

F

H

H

bleue

noir

gris

bleue

bleue

gris

bleue

noir

noir

gris

DataFrame

df=data.frame(T, P, S, C)
head(df)

II- Analyse descriptive univaries

---

Lecture des données

#Lecture des données dans le fichier DBM1.csv
data=read.csv("Ressources/DBM1.csv", sep = ";")

#Ou utiliser la commande et ne plus convertir en facteur les variables categorielles
data=read.csv("Ressources/DBM1.csv", sep=";", stringsAsFactors = TRUE)

#Conversion en facteur des variables categorielles
data$S=as.factor(data$S)
data$C=as.factor(data$C)

head(data)

Parametres de position

1- Moyenne

mean(data$T)

[1] 174.3111

mean(data$P)

[1] 65.08889

2- Mediane

median(data$T)

[1] 177

median(data$F)

[1] 2

3- Quantiles

quantile(data$T)

  0%  25%  50%  75% 100% 
 157  167  177  181  190 

quantile(data$F)

  0%  25%  50%  75% 100% 
   0    1    2    3    7 

quantile(data$P, probs = 0.8)

80% 
 72 

4- Valeurs extremes (min, max)

min(data$F)

[1] 0

max(data$P)

[1] 98

5- Mode

library(conflicted)
#install.packages("DescTools")
library(DescTools)

Warning: le package 'DescTools' a été compilé avec la version R 4.3.3

Mode(data$S, na.rm = TRUE)

[1] h
attr(,"freq")
[1] 30
Levels: f h

Mode(data$C)

[1] brun
attr(,"freq")
[1] 23
Levels: bleu brun gris noir vert

6- Resume descriptif

summary(data)

       T               P         S            F            C     
 Min.   :157.0   Min.   :47.00   f:15   Min.   :0.000   bleu:12  
 1st Qu.:167.0   1st Qu.:59.00   h:30   1st Qu.:1.000   brun:23  
 Median :177.0   Median :65.00          Median :2.000   gris: 1  
 Mean   :174.3   Mean   :65.09          Mean   :1.933   noir: 2  
 3rd Qu.:181.0   3rd Qu.:71.00          3rd Qu.:3.000   vert: 7  
 Max.   :190.0   Max.   :98.00          Max.   :7.000            

#install.packages("psych")
library(psych)

Warning: le package 'psych' a été compilé avec la version R 4.3.3

describe(data[,c("T", "P", "F")])

Desc(data)

------------------------------------------------------------------------------ 
Describe data (data.frame):

data frame: 45 obs. of  5 variables
        45 complete cases (100.0%)

  Nr  ColName  Class    NAs  Levels                                     
  1   T        integer  .                                               
  2   P        integer  .                                               
  3   S        factor   .    (2): 1-f, 2-h                              
  4   F        integer  .                                               
  5   C        factor   .    (5): 1-bleu, 2-brun, 3-gris, 4-noir, 5-vert


------------------------------------------------------------------------------ 
1 - T (integer)

  length       n     NAs  unique      0s    mean  meanCI'
      45      45       0      26       0  174.31  171.66
          100.0%    0.0%            0.0%          176.97
                                                        
     .05     .10     .25  median     .75     .90     .95
  160.20  162.40  167.00  177.00  181.00  184.00  185.00
                                                        
   range      sd   vcoef     mad     IQR    skew    kurt
   33.00    8.84    0.05   10.38   14.00   -0.26   -1.20
                                                        
lowest : 157, 158, 160, 161, 162
highest: 183 (4), 184 (2), 185 (2), 188, 190

heap(?): remarkable frequency (15.6%) for the mode(s) (= 180)

' 95%-CI (classic)

------------------------------------------------------------------------------ 
2 - P (integer)

  length       n    NAs  unique     0s   mean  meanCI'
      45      45      0      21      0  65.09   62.40
          100.0%   0.0%           0.0%          67.78
                                                     
     .05     .10    .25  median    .75    .90     .95
   52.00   54.60  59.00   65.00  71.00  74.20   75.00
                                                     
   range      sd  vcoef     mad    IQR   skew    kurt
   51.00    8.95   0.14    8.90  12.00   0.74    2.27
                                                     
lowest : 47, 49, 52 (2), 53, 57 (3)
highest: 72 (5), 73, 75 (3), 78, 98

heap(?): remarkable frequency (13.3%) for the mode(s) (= 65)

' 95%-CI (classic)

------------------------------------------------------------------------------ 
3 - S (factor - dichotomous)

  length      n    NAs unique
      45     45      0      2
         100.0%   0.0%       

   freq   perc  lci.95  uci.95'
f    15  33.3%   21.4%   47.9%
h    30  66.7%   52.1%   78.6%

' 95%-CI (Wilson)

------------------------------------------------------------------------------ 
4 - F (integer)

  length       n    NAs  unique     0s  mean  meanCI'
      45      45      0       8      8  1.93    1.45
          100.0%   0.0%          17.8%          2.41
                                                    
     .05     .10    .25  median    .75   .90     .95
    0.00    0.00   1.00    2.00   3.00  4.00    4.80
                                                    
   range      sd  vcoef     mad    IQR  skew    kurt
    7.00    1.60   0.83    1.48   2.00  1.05    1.08
                                                    

   value  freq   perc  cumfreq  cumperc
1      0     8  17.8%        8    17.8%
2      1    12  26.7%       20    44.4%
3      2    12  26.7%       32    71.1%
4      3     7  15.6%       39    86.7%
5      4     3   6.7%       42    93.3%
6      5     1   2.2%       43    95.6%
7      6     1   2.2%       44    97.8%
8      7     1   2.2%       45   100.0%

' 95%-CI (classic)

------------------------------------------------------------------------------ 
5 - C (factor)

  length      n    NAs unique levels  dupes
      45     45      0      5      5      y
         100.0%   0.0%                     

   level  freq   perc  cumfreq  cumperc
1   brun    23  51.1%       23    51.1%
2   bleu    12  26.7%       35    77.8%
3   vert     7  15.6%       42    93.3%
4   noir     2   4.4%       44    97.8%
5   gris     1   2.2%       45   100.0%

Parametre de dispersion

1- Variance

var(data$T)

[1] 78.08283

var(data$P)

[1] 80.17374

var(data$F)

[1] 2.563636

2- Ecart Type

sd(data$T)

[1] 8.836449

sd(data$P)

[1] 8.953979

sd(data$F)

[1] 1.601136

3- Etendue

max(data$T)-min(data$T)

[1] 33

max(data$P)-min(data$P)

[1] 51

max(data$F)-min(data$F)

[1] 7

Visualisation des données

#install.packages("ggplot2")
library(ggplot2)

Warning: le package 'ggplot2' a été compilé avec la version R 4.3.3

1- Variables Numeriques

ggplot(data, aes(x=F))+geom_histogram(bins = length(unique(data$F)), col="black",  fill = "lightgray")+
  labs(x = "Nombre de frere", y = "Frequence", title="Distribution suivant le nombre de freres")

ggplot(data, aes(x=F))+geom_boxplot( fill = "lightgray")+
  labs(x = "Nombre de freres", title="Boxplot du nombre de freres des individus")

ggplot(data, aes(x=P))+geom_boxplot( fill = "lightgray")+
  labs(x = "Poids", title="Boxplot du poids des individus")

ggplot(data, aes(x=T))+geom_boxplot( fill = "lightgray")+
  labs(x = "Taille", title="Boxplot de la taille des individus")

ggplot(data, aes(x=T))+geom_histogram(aes(y=..density..),fill="lightgray", col="black")+
  geom_density(linewidth=0.55, col="black")+
  labs(x = "Taille", y = "Frequence", title="Distribution suivant la taille")

Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(density)` instead.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
generated.

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

2- Variables Categorielles

ggplot(data, aes(x=S))+geom_bar(fill="chocolate")+
  labs(x = "Sexe", y = "Frequence", title="Distribution suivant le sexe")

tab_C=table(data$C)

#On recupere les lables 
labs_C=rownames(tab_C)

labs_C=paste(labs_C, round(100*tab_C/sum(tab_C), 2), sep=" (")
labs_C=paste(labs_C, "", sep="%)")
pie(tab_C, 
    col=c("#999999", "#E69F00", "#56B4E9","#E9967A","#FFA07A"),
    main="Distribution de la couleur des yeux")

legend("topright", labs_C, cex=0.8,fill=rainbow(length(tab_C)))

#install.packages("plotrix")
library(plotrix)

pie3D(tab_C, explode=0.08,
      border = "white",
      height = 0.25
      ,main="Distribution de la couleur des yeux")
legend("topright", labs_C, cex=0.5,fill=rainbow(length(tab_C)), 
       ncol = 1)

III- Analyses descriptives bivaries

---

1- Quanti-Quanti

#Covariance avec la methode de pearson
print("pearson")

[1] "pearson"

cov(data[,c("T","P","F")], method ="pearson" )

          T         P         F
T 78.082828 45.698990 -2.001515
P 45.698990 80.173737 -4.334848
F -2.001515 -4.334848  2.563636

#Covariance avec la methode de kendall
print("kendall")

[1] "kendall"

cov(data[,c("T","P","F")], method ="kendall" )

     T    P    F
T 1902  764 -214
P  764 1888 -392
F -214 -392 1612

#Covariance avec la methode de spearman
print("spearman")

[1] "spearman"

cov(data[,c("T","P","F")], method ="spearman" )

          T        P         F
T 171.59091 104.5227 -27.46591
P 104.52273 171.5568 -51.75000
F -27.46591 -51.7500 164.36364

#Correlation avec la methode de pearson
print("pearson")

[1] "pearson"

cor(data[,c("T","P","F")], method ="pearson" )

           T          P          F
T  1.0000000  0.5775808 -0.1414663
P  0.5775808  1.0000000 -0.3023637
F -0.1414663 -0.3023637  1.0000000

#Correlation avec la methode de kendall
print("kendall")

[1] "kendall"

cor(data[,c("T","P","F")], method ="kendall" )

           T          P          F
T  1.0000000  0.4031690 -0.1222154
P  0.4031690  1.0000000 -0.2246997
F -0.1222154 -0.2246997  1.0000000

#Correlation avec la methode de spearman
print("spearman")

[1] "spearman"

cor(data[,c("T","P","F")], method ="spearman" )

           T          P          F
T  1.0000000  0.6091996 -0.1635475
P  0.6091996  1.0000000 -0.3081793
F -0.1635475 -0.3081793  1.0000000

library(heatmaply)

Le chargement a nécessité le package : plotly

Warning: le package 'plotly' a été compilé avec la version R 4.3.3

Attachement du package : 'plotly'

L'objet suivant est masqué depuis 'package:ggplot2':

    last_plot

L'objet suivant est masqué depuis 'package:stats':

    filter

L'objet suivant est masqué depuis 'package:graphics':

    layout

Le chargement a nécessité le package : viridis

Warning: le package 'viridis' a été compilé avec la version R 4.3.3

Le chargement a nécessité le package : viridisLite

Registered S3 methods overwritten by 'registry':
  method               from 
  print.registry_field proxy
  print.registry_entry proxy

======================
Welcome to heatmaply version 1.5.0

Type citation('heatmaply') for how to cite the package.
Type ?heatmaply for the main documentation.

The github page is: https://github.com/talgalili/heatmaply/
Please submit your suggestions and bug-reports at: https://github.com/talgalili/heatmaply/issues
You may ask questions at stackoverflow, use the r and heatmaply tags: 
     https://stackoverflow.com/questions/tagged/heatmaply
======================

#Visualisation de la matrice de correlation
heatmaply(cor(data[,c("T","P","F")], method ="spearman" ), 
          main = "Matrice de correlation")

#Visualisation de la matrice de covariance
heatmaply(cov(data[,c("T","P","F")], method ="spearman" ), 
          main = "Matrice de covariance")

#Nuage de points
model=lm(T~P, data)
coefs=coef(model)
ggplot(data, aes(x=P, y=T))+geom_point(col="blue")+
  geom_abline(intercept = coefs[1], slope=coefs[2], linewidth=1)+
  labs(title =)+
  labs(title="Taille en fonction du poids", x="Taille",
      y="Poids")

#Nuage de points
ggplot(data, aes(x=P, y=F))+geom_point(col="blue")+
  geom_smooth(method = "lm")+
  labs(title="Nombre de frere en fonction du poids", x="Poids",
      y="Nombre de frere")

`geom_smooth()` using formula = 'y ~ x'

#Nuage de points
ggplot(data, aes(x=T, y=F))+geom_point(col="blue")+
  geom_smooth(method = "lm")+
  labs(title="Nombre de frere en fonction de la taille", x="Taille",
      y="Nombre de frere")

`geom_smooth()` using formula = 'y ~ x'

#install.packages("GGally")
library(GGally)

Warning: le package 'GGally' a été compilé avec la version R 4.3.3

Registered S3 method overwritten by 'GGally':
  method from   
  +.gg   ggplot2

ggcorr(data[,c("P","T","F")])

2-Quali-Quali

table(data$S,data$C)

   
    bleu brun gris noir vert
  f    4    8    0    0    3
  h    8   15    1    2    4

ggplot(data, aes(x=C))+geom_bar(aes(fill=S), position = "dodge")+
  scale_fill_brewer(palette = "Pastel2")+
  labs(title = "Distribution des individus suivant le sexe et la couleur des yeux",
       x="Couleur des yeux", y="Frequence",fill="sexe")

library(grid)
library(vcd)

Warning: le package 'vcd' a été compilé avec la version R 4.3.3

#Diagramme en mosaic
mosaicplot(~ C + S , 
           data = data, 
           color = TRUE,
           las=1,
           main = "Diagramme en mosaic du sexe et de la couleur des yeux")

#Matrice de contengence
ggplot(data, aes(y=S,x=C))+geom_bin_2d()+
  labs(title = "Distribution des individus suivant le sexe et la couleur des yeux",
       x="Couleur des yeux", y="sexe")

#Parametres sur les variables qualitatives
ass=assocstats(table(data$S, data$C))
print(ifelse(ass$cramer<=0.3,
             paste("v cramer:",round(ass$cramer,2),"Le sexe et faiblement associé à la couleur des yeux", sep=" "),
             ifelse(ass$cramer>0.3 & ass$cramer<=0.6,paste("v cramer:",round(ass$cramer,2),"Le sexe et moyennement associé à la couleur des yeux", sep = " "),
              paste("v cramer:",round(ass$cramer,2),"Le sexe et moyennement associé à la couleur des yeux", sep = " "))))

[1] "v cramer: 0.2 Le sexe et faiblement associé à la couleur des yeux"

3- Quali-Quanti

ggplot(data, aes(x=P, y=S))+geom_violin(aes(fill=S))+
  labs(x = "Poids", y="Sexe", fill="Sexe",
       title="Violinplot du poids des individus suivant le sexe")+
  scale_fill_brewer(palette = "Pastel1")

ggplot(data, aes(y=P, x=C))+geom_boxplot(aes(fill=C),outlier.colour = "red",)+
  labs(y = "Poids", x="Couleur des yeux", fill="Couleur des yeux",
       title="Boxplot de la taille des individus suivant la couleur des yeux")+
  scale_fill_brewer(palette = "Pastel2")

ggpairs(data[,c("P","T","F")])

ggally_density(data[,c("P","T","S")],aes(x=P, y=T, col=S))+labs(x="Poids", y="Taille", col="Sexe")

ggplot(data, aes(P, T))+stat_density_2d(aes(fill=S), geom = "polygon")+labs(x="Poids", y="Taille", col="Sexe")

IV- Statistiques Inférentielles

---

On rejete H0 losque p-value < α (5% par defaut)

Test d’Homogénéité des Variances

Pour tester l’homogeneite des variances (égalité des variances), de nombreux tests peuvent être utilisés notament :

• Test F : Comparez les variances de deux groupes. Les données doivent être normalement distribuées.

• Test de Bartlett : Comparer les variances de deux groupes ou plus. Les données doivent être normalement distribuées.

• Le test de Levene : Une alternative robuste au test de Bartlett qui est moins sensible aux écarts de normalité.

• Test de Fligner-Killeen : un test non paramétrique qui est très robuste contre les écarts de normalité.

NB : Il est à noter que le test de Levene est le plus couramment utilisé dans la littérature.

1- Test F

#install.packages("tidyverse")
#install.packages("conflicted")
#install.packages("rstatix")
library(conflicted)
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ lubridate 1.9.3     ✔ tibble    3.2.1
✔ purrr     1.0.2     ✔ tidyr     1.3.0

library(rstatix)

Warning: le package 'rstatix' a été compilé avec la version R 4.3.3

library(ggpubr)

#On visualise la normalité de la variable
ggqqplot(data, x="P")+labs(title="qqplot du poids de l'ensemble des individus")

ggqqplot(data, x="P", facet.by = "S")+labs(title="qqplot du poids de chaque sexe")

Les points étant presque tous situés autour de la droite de reference, nous pouvons supposer une normalité pour la variable poids. Cette normalité est vérifié par le test de shapiro wilk.

data %>% group_by(S) %>% shapiro_test(P)

NB : Si la taille de votre échantillon est supérieure à 50, le tracé QQ normal est préféré car avec des échantillons plus grands, le test de Shapiro-Wilk devient très sensible même à un écart mineur par rapport à la normalité.

var.test(P~S, data)

    F test to compare two variances

data:  P by S
F = 0.53037, num df = 14, denom df = 29, p-value = 0.2111
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.2251757 1.4533624
sample estimates:
ratio of variances 
         0.5303741 

La P-Value > 5% indique que l’on accepte l’hypothèse H0 et donc que les variances des poids des deux sexes (h et f) sont égales.

2- Comparer plusieurs écarts

Cette section décrit comment comparer plusieurs variances dans R à l’aide des tests de Bartlett, Levene ou Fligner-Killeen.

Hypothèses statistiques . Pour tous ces tests qui suivent, l’hypothèse nulle est que les variances de toutes les populations sont égales, l’hypothèse alternative est qu’au moins deux d’entre elles diffèrent. Par conséquent, des valeurs p inférieures à 0,05 suggèrent que les variances sont significativement différentes et que l’hypothèse d’homogénéité de la variance a été violée.

2.1- test de Bartlett

bartlett.test(P~S, data=data)

    Bartlett test of homogeneity of variances

data:  P by S
Bartlett's K-squared = 1.6963, df = 1, p-value = 0.1928

2.2- Test de Levene

LeveneTest(P~S, data=data)

2.3- Test de Fligner-Killeen

fligner.test(P~S, data=data)

    Fligner-Killeen test of homogeneity of variances

data:  P by S
Fligner-Killeen:med chi-squared = 0.20196, df = 1, p-value = 0.6531

T test

Ce test permet de comparer les moyennes de deux groupes ou d’un groupe avec une valeur de reference connue

1- t test à un échantillon

Dans cette situation nous avons un echantillon {x1, x2, …., xn} et on veut tester si la moyenne de cette echantillons est egale, superieur ou inferieur a une valeur connue.

On affirme que le poids moyen de la population est de 65kg. On utilise un t test pour verifier cette hypothese.

H0: mu=65

H1: mu<>65

mean(data$P)

[1] 65.08889

res<-t.test(data$P, mu = 65)
res

    One Sample t-test

data:  data$P
t = 0.066594, df = 44, p-value = 0.9472
alternative hypothesis: true mean is not equal to 65
95 percent confidence interval:
 62.39882 67.77896
sample estimates:
mean of x 
 65.08889 

On observe que la p-value > 5% donc on a bien le poids moyen de la population égale à 65 avec 5% de chance de se tromper.

res<-t_test(data, P~1, mu=65)
# Create the box-plot
bxp <- ggboxplot(
  data$P, width = 0.5, add = c("mean", "jitter"), 
  ylab = "Weight (g)", xlab = FALSE
  )

Warning: The `fun.y` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
ℹ Please use the `fun` argument instead.
ℹ The deprecated feature was likely used in the ggpubr package.
  Please report the issue at <https://github.com/kassambara/ggpubr/issues>.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
generated.

Warning: The `fun.ymin` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
ℹ Please use the `fun.min` argument instead.
ℹ The deprecated feature was likely used in the ggpubr package.
  Please report the issue at <https://github.com/kassambara/ggpubr/issues>.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
generated.

Warning: The `fun.ymax` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
ℹ Please use the `fun.max` argument instead.
ℹ The deprecated feature was likely used in the ggpubr package.
  Please report the issue at <https://github.com/kassambara/ggpubr/issues>.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
generated.

# Add significance levels
bxp + labs(subtitle = get_test_label(res, detailed = TRUE))

De même on peut tester si la taille moyenne de la population est superieure a 150cm

H0: mu=150

H1: mu > 150

mean(data$T)

[1] 174.3111

t.test(data$T, mu=150, alternative = "greater")

    One Sample t-test

data:  data$T
t = 18.456, df = 44, p-value < 2.2e-16
alternative hypothesis: true mean is greater than 150
95 percent confidence interval:
 172.0978      Inf
sample estimates:
mean of x 
 174.3111 

On observe que la p-value < 5%. Donc avec un risque de 5% de se tromper la taille de la population est superieur a 150cm

En utilisant cette fois le package rstatix on peut aussi tester pareillement si le nombre moyen de frere est inferieur a 3

H0: mu=3

H1: mu < 3

stat.test=data %>% t_test(F~1, mu=3, alternative = "less", detailed = TRUE)
stat.test

On observe ici qu’avec une confiance de 95% que le nombre de frère est bien inferieur à 3.

ggdensity(data, x = "F", rug = TRUE, fill = "lightgray") +
  stat_central_tendency(type = "mean", color = "red", linetype = "dashed") +
  geom_vline(xintercept = 3 , color = "blue", linetype = "dashed") + 
  labs(subtitle = get_test_label(stat.test,  detailed = TRUE), x="Nombre de frere")

2- t test pour un échantillon indépendant

Le test t pour les écahntillons indépendants permet comme son nom l’indique de comparer les moyennes de deux échantillons biens différents. Il est aussi appelé t test pour un ecahntillon apparié.

Prérequis :

• Une variable dependante catégorielle

• Une variable indépendante numérique

• Indépendance entre les observations de chaque catégorie

• Données normalement distribuées

• Vérifier l’égalité des variances (en fonction des résultats du test de levene mettre var.equal = TRUE/FALSE)

data %>% group_by(S) %>% 
  get_summary_stats(P, type = "mean_sd")

On peut supposer que le poids des hommes est superieur a celui des femme

H0: mu(h)=mu(f)

H1: mu > mu(f)

c’est a dire tester si

H0: mu(h)-mu(f)=0

H1: mu - mu(f)>0

#var.equal = TRUE car les variances sont égales d'après le test de levene
t.test(P~S, data=data, alternative="less",
       var.equal = TRUE)

    Two Sample t-test

data:  P by S
t = -4.4418, df = 43, p-value = 3.074e-05
alternative hypothesis: true difference in means between group f and group h is less than 0
95 percent confidence interval:
      -Inf -6.546832
sample estimates:
mean in group f mean in group h 
       58.06667        68.60000 

Donc on peut affirmer avec 95% de chance qu’en moyenne les hommes pesent plus que les femmes

res<-t_test(data, P~S, alternative="less", var.equal = TRUE)
res <- res %>% add_xy_position(x = "S")
# Create the box-plot
bxp <- ggboxplot(
  data, x="S", y="P", width = 0.5, add = c("mean"), 
  ylab = "Poids (kg)", xlab = "Sexe"
  )
# Add significance levels
bxp + 
  stat_pvalue_manual(res, tip.length = 0) +
  labs(subtitle = get_test_label(res, detailed = TRUE), title = "Boxplot du poids en fonction du sexe")

On peut évaluer cette différence en utilisant le d de cohens.

3- t test par pair

Ici on veut comparer les moyennes de deux mesures effectués sur des mêmes individus. Par exemple lorsqu’on veut comparer les mesures de tension arterielles des patients avant et après un traitement X.

#install.packages("datarium")
library(datarium)

Warning: le package 'datarium' a été compilé avec la version R 4.3.3

data("mice2", package = "datarium")
head(mice2, 3)

#Transformation des données en une colonne pour la classe et une pour les données mesurées
mice2.long <- mice2 %>%
  gather(key = "group", value = "weight", before, after)
head(mice2.long)

#Resumé des données 
mice2.long %>%
  group_by(group) %>%
  get_summary_stats(weight, type = "mean_sd")

#Execution du test
res<-t_test(weight~group, data=mice2.long, paired = TRUE, detailed = TRUE)
res

On peut conclure que le poids est significativement différent avant et après le traitement car la p-value est < 5%.

On peut évaluer cette différence en utilisant le d de cohens comme suit :

cohens_d(weight~group, data=mice2.long, paired = TRUE)

Donc l’écart entres les poids moyennes avant et après est d’environ 8.08

bxp <- ggpaired(mice2.long, x = "group", y = "weight", 
         order = c("before", "after"),
         ylab = "Weight", xlab = "Groups")

res <- res %>% add_xy_position(x = "group")
bxp + 
  stat_pvalue_manual(res, tip.length = 0) +
  labs(subtitle = get_test_label(res, detailed= TRUE))

Test binomial

Ce test d’applique a des variables binomiales telles que le sexe (homme et femme). Il permet de comparer une proportion d’un échantillon (exemple la proportion des hommes) avec une valeur de reference.

Exemple :

On suppose qu’il ya autant d’homme que de femme dans la population

H0 : P(h)=P(f)

H1 : P(h)<>P(f)

table(data$S)

 f  h 
15 30 

binom.test(table(data$S))

    Exact binomial test

data:  table(data$S)
number of successes = 15, number of trials = 45, p-value = 0.0357
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
 0.2000057 0.4895036
sample estimates:
probability of success 
             0.3333333 

p-value < 5% on rejete l’hypothèse nulle et donc la proportion des hommes est significatvement différente de celle des femmes

Test de khi 2

Le test du Khi-deux est un test d’hypothèse utilisé pour déterminer s’il existe une relation entre deux variables catégorielles.

data$PR<-cut(data$P, breaks=5)
ggplot(data, aes(x=PR, fill=S))+geom_bar(position = "dodge")

On pourrait supposer une liaison entre le sexe et la tranche de poids

khi_test<-chisq.test(x=data$S, y=data$PR)

Warning in chisq.test(x = data$S, y = data$PR): L’approximation du Chi-2 est
peut-être incorrecte

print("Distribution observe")

[1] "Distribution observe"

khi_test$observed

      data$PR
data$S (46.9,57.2] (57.2,67.4] (67.4,77.6] (77.6,87.8] (87.8,98.1]
     f           6           8           1           0           0
     h           2          11          15           1           1

print("Distribution attendue")

[1] "Distribution attendue"

khi_test$expected

      data$PR
data$S (46.9,57.2] (57.2,67.4] (67.4,77.6] (77.6,87.8] (87.8,98.1]
     f    2.666667    6.333333    5.333333   0.3333333   0.3333333
     h    5.333333   12.666667   10.666667   0.6666667   0.6666667

khi_test

    Pearson's Chi-squared test

data:  data$S and data$PR
X-squared = 13.189, df = 4, p-value = 0.01039

NB : Les distributions théoriques (attendue) doivent toutes être supérieures à 5

Test U de Mann-Whitney

Le test U de Mann-Whitney peut être utilisé pour vérifier s’il existe une différence entre deux échantillons (groupes), et les données ne doivent pas nécessairement être distribuées normalement. Pour pouvoir calculer un test U de Mann-Whitney, il faut disposer de deux échantillons aléatoires indépendants présentant au moins des caractéristiques à échelle ordinale.

H0 : Il n’existe pas de différence entre les groupes (du point de vu des tendances centrales)

group_by(data,S) %>%
  summarise(
    count = n(),
    median = median(P, na.rm = TRUE),
    IQR = IQR(P, na.rm = TRUE))

wilcox.test(P~S, data = data)

Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): impossible
de calculer la p-value exacte avec des ex-aequos

    Wilcoxon rank sum test with continuity correction

data:  P by S
W = 56, p-value = 4.739e-05
alternative hypothesis: true location shift is not equal to 0

Test de Wilcoxon

Le test de Wilcoxon (test du rang de signe de Wilcoxon) vérifie si les valeurs moyennes de deux groupes dépendants diffèrent significativement l’une de l’autre.

Le test de Wilcoxon est un test non paramétrique et est donc soumis à beaucoup moins de conditions prélables que son homologue paramétrique, le test t pour échantillons dépendants. Par conséquent, dès que les conditions limites du test t pour échantillons dépendants ne sont plus remplies, le test de Wilcoxon est utilisé.

Conditions du test de Wilcoxon

Le test de Wilcoxon étant un test non paramétrique, il n’est pas nécessaire que les données soient normalement distribuées. Toutefois, pour calculer un test de Wilcoxon, les échantillons doivent être dépendants. Les échantillons dépendants sont présents, par exemple, lorsque les données sont obtenues à partir de mesures répétées ou lorsqu’il s’agit de paires dites naturelles.

• Mesure répétée : une caractéristique d’une personne, par exemple son poids, a été mesurée à deux moments différents.

• Couples naturels : les valeurs ne doivent pas nécessairement provenir de la même personne, mais de personnes qui vont ensemble, par exemple avocat/client, épouse/mari et psychologue/patient. Bien entendu, il ne s’agit pas nécessairement de personnes.

• Indépendance : le test de Wilcoxon suppose l’indépendance, c’est-à-dire que les observations appariées sont tirées au hasard et de manière indépendante.

En outre, la forme de la distribution des différences entre les deux échantillons dépendants doit être approximativement symétrique.

Si les données ne sont pas disponibles par paires, le test U de Mann-Whitney est utilisé à la place du test de Wilcoxon.

# Load the data from datarium package
data("genderweight", package = "datarium")
# Show a sample of the data by group
set.seed(123)
genderweight %>% sample_n_by(group, size = 2)

genderweight %>%
  group_by(group) %>%
  get_summary_stats(weight, type = "median_iqr")

bxp <- ggboxplot(
  genderweight, x = "group", y = "weight", 
  ylab = "Weight", xlab = "Groups", add = c("jitter") 
  )
bxp

stat.test <- genderweight %>% 
  wilcox_test(weight ~ group) %>%
  add_significance()
stat.test

stat.test <- stat.test %>% add_xy_position(x = "group")
bxp + 
  stat_pvalue_manual(stat.test, tip.length = 0) +
  labs(subtitle = get_test_label(stat.test, detailed = TRUE))

Anova

Le test ANOVA (ou Analyse de Variance en francais) est utilisé pour comparer la moyenne de plusieurs groupes. Le terme ANOVA est un peu trompeur. Bien que le nom de la technique fasse référence aux variances, l’objectif principal de l’ANOVA est d’étudier les différences de moyennes.

1- Anova a un facteur

Une extension du test t pour échantillons indépendants pour comparer les moyennes dans une situation où il y a plus de deux groupes. Il s’agit du cas le plus simple de test ANOVA où les données sont organisées en plusieurs groupes selon une seule variable de regroupement (également appelée variable factorielle). D’autres synonymes sont : ANOVA à 1 voie , ANOVA à un facteur et ANOVA inter-sujets 

data("PlantGrowth")
PlantGrowth <- PlantGrowth %>%
  reorder_levels(group, order = c("ctrl", "trt1", "trt2"))

PlantGrowth %>%
  group_by(group) %>%
  get_summary_stats(weight, type = "mean_sd")

ggboxplot(PlantGrowth, x = "group", y = "weight")

Les valeurs aberrantes peuvent être facilement identifiées à l’aide des méthodes de diagramme en boîte, implémentées dans la fonction R identify_outliers()du package rstatix.

PlantGrowth %>% 
  group_by(group) %>%
  identify_outliers(weight)

Vérification de l’hypothèse de normalité.
Elle peut être fait soit :

• A partir des tests de normalité dans chaque groupe

• Ou a partir des résidus du modèle ANOVA

# Construction du modèle linéaire
model  <- lm(weight ~ group, data = PlantGrowth)
# QQ Plot
ggqqplot(residuals(model))

Les valeurs sont proches de la droite de reférence, cela suppose une normalité des résidus

H0 : La variable est normalement distribué

#Test de normalité de shapiro wilk 
shapiro_test(residuals(model))

Dans le cas de la normalité de chaque groupe on aurait

ggqqplot(data=PlantGrowth, x="weight", facet.by = "group")

PlantGrowth %>% group_by(group) %>%
  shapiro_test(weight)

La P-value ici est > 5% donc le poids suit une distribution normale.

res.aov <-PlantGrowth %>% anova_test(weight ~ group)

2- Tests post-hoc

Une ANOVA unidirectionnelle significative est généralement suivie de tests post-hoc de Tukey pour effectuer plusieurs comparaisons par paires entre les groupes. Fonction de la touche R : tukey_hsd()[rstatix].

pwc <- PlantGrowth %>% tukey_hsd(weight ~ group)

#Visualisation 
pwc <- pwc %>% add_xy_position(x = "group")
ggboxplot(PlantGrowth, x = "group", y = "weight") +
  stat_pvalue_manual(pwc, hide.ns = TRUE) +
  labs(
    subtitle = get_test_label(res.aov, detailed = TRUE),
    caption = get_pwc_label(pwc)
    )

NB : Le test ANOVA unidirectionnel classique nécessite une hypothèse de variances égales pour tous les groupes. Dans le cas contraire, il est conseillé d’utiliser :

• Le test unidirectionnel de Welch est une alternative à l’ANOVA unidirectionnelle standard dans les situations où l’homogénéité de la variance ne peut pas être supposée (c’est-à-dire que le test de Levene est significatif).

• Dans ce cas, le test post hoc de Games-Howell ou les tests t par paires (sans hypothèse de variances égales) peuvent être utilisés pour comparer toutes les combinaisons possibles de différences de groupe.

#Dans ce cas voici les fonctions à utiliser
#welch_anova_test()
#games_howell_test()

3- Anova à deux facteurs

Utilisée pour évaluer simultanément l’effet de deux variables de regroupement différentes sur une variable de résultat continue. D’autres synonymes sont : deux plans factoriels , anova factorielle ou ANOVA bidirectionnelle entre sujets .

data("jobsatisfaction", package = "datarium")
jobsatisfaction %>%
  group_by(gender, education_level) %>%
  get_summary_stats(score, type = "mean_sd")

ggplot(jobsatisfaction)+geom_boxplot(aes(x=gender, y=score, fill=education_level))

model  <- lm(score ~ gender*education_level,
             data = jobsatisfaction)
ggqqplot(residuals(model))

shapiro_test(residuals(model))

jobsatisfaction %>%
  group_by(gender, education_level) %>%
  shapiro_test(score)

ggqqplot(jobsatisfaction, "score", ggtheme = theme_bw()) +
  facet_grid(gender ~ education_level)

jobsatisfaction %>% levene_test(score ~ gender*education_level)

res.aov <- jobsatisfaction %>% anova_test(score ~ gender * education_level)
res.aov

4- Tests post-hoc

Une interaction bidirectionnelle significative indique que l’impact d’un facteur (par exemple, le niveau_d’éducation) sur la variable de résultat (par exemple, le score de satisfaction au travail) dépend du niveau de l’autre facteur (par exemple, le sexe) (et vice versa). Ainsi, vous pouvez décomposer une interaction bidirectionnelle significative en :

• Effet principal simple : exécuter un modèle unidirectionnel de la première variable à chaque niveau de la deuxième variable,

• Comparaisons simples par paires : si l’effet principal simple est significatif, effectuez plusieurs comparaisons par paires pour déterminer quels groupes sont différents.

Pour une interaction bidirectionnelle non significative , vous devez déterminer si vous avez des effets principaux statistiquement significatifs à partir du résultat de l’ANOVA. Un effet principal significatif peut être suivi par des comparaisons par paires entre les groupes.

model <- lm(score ~ gender * education_level, data = jobsatisfaction)
jobsatisfaction %>%
  group_by(gender) %>%
  anova_test(score ~ education_level, error = model)

L’effet principal simple du « niveau_d’éducation » sur le score de satisfaction au travail était statistiquement significatif tant pour les hommes que pour les femmes (p < 0,0001).

En d’autres termes, il existe une différence statistiquement significative dans le score moyen de satisfaction au travail entre les hommes ayant fait des études scolaires, collégiales ou universitaires, F(2, 52) = 132, p < 0,0001. La même conclusion est vraie pour les femmes , F(2, 52) = 62,8, p < 0,0001.

NB : Il convient de noter que la signification statistique des analyses simples de l’effet principal a été acceptée à un niveau alpha ajusté par Bonferroni de 0,025. Cela correspond au niveau actuel auquel vous déclarez la signification statistique (c’est-à-dire p < 0,05) divisé par le nombre d’effets principaux simples que vous calculez (c’est-à-dire 2).

V- Analyse multivaries

---

ACP : Analyse en Composantes Principales

AFC : Analyse Factorielle de Correspondance

ACM : Analyse Factorielle de Correspondance Multiple

CAH : Classification Ascendante Hierarchique

VI- Documentation

---

https://datatab.fr/tutorial/get-started

https://3mmarand.github.io/comp4biosci/

https://epirhandbook.com/en/index.html

https://openintro-ims.netlify.app/

Realiser l’ANOVA avec R

https://larmarange.github.io/analyse-R/graphiques-bivaries-ggplot2.html

https://www.datanovia.com/en/lessons/how-to-do-a-t-test-in-r-calculation-and-reporting/

https://www.umr-lastig.fr/paul-chapron/resources/cours_site/galerie-de-graphiques-avec-ggplot.html

https://www.datanovia.com/en/fr/lessons/test-dhomogeneite-des-variances-dans-r/

https://jokergoo.github.io/ComplexHeatmap-reference/book/introduction.html

http://www.sthda.com/english/wiki/descriptive-statistics-and-graphics

http://www.sthda.com/english/wiki/ggpubr-create-easily-publication-ready-plots